home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_6.1 / xvor / geometry.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  4.9 KB  |  246 lines
  1. /* 
  2.  * geometry.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * GEOMETRY.C
  11.  *
  12.  */
  13. #include <stdio.h>
  14. #include <stdlib.h>
  15. #include <math.h>
  16. #include "vor.h"
  17.  
  18. #define    ZERO    0.0
  19.  
  20.  
  21. void
  22. geominit ()
  23. {
  24.   float sn;
  25.  
  26.   freeinit (&efl, sizeof (struct Edge));
  27.  
  28.   nedges = nvertices = 0;
  29.   sn = (float) (nsites + 4);    /* significance of const. 4 ?? */
  30.   sqrt_nsites = 1 + (int) sqrt (sn);  /* sqrt returns double !! */
  31.  
  32.   deltay = ymax - ymin;
  33.   deltax = xmax - xmin;
  34. }
  35.  
  36.  
  37. struct Edge *
  38. bisect (s1, s2, imgIO, value)
  39.      struct Site *s1, *s2;
  40.      Image *imgIO;
  41.      int value;
  42. {
  43.   double dx, dy, adx, ady;
  44.   struct Edge *newedge;
  45.  
  46.   newedge = (struct Edge *) getfree (&efl);
  47.  
  48.   newedge->reg[0] = s1;
  49.   newedge->reg[1] = s2;
  50.   ref (s1);
  51.   ref (s2);
  52.   newedge->ep[0] = (struct Site *) NULL;
  53.   newedge->ep[1] = (struct Site *) NULL;
  54.  
  55.   dx = s2->coord.x - s1->coord.x;
  56.   dy = s2->coord.y - s1->coord.y;
  57.   if ((-1.0e-10 < dx) && (dx < 1.0e-10))
  58.     dx = ZERO;                  /* vert */
  59.   if ((-1.0e-10 < dy) && (dy < 1.0e-10))
  60.     dy = ZERO;                  /* hor */
  61.   adx = dx > 0 ? dx : -dx;
  62.   ady = dy > 0 ? dy : -dy;
  63.   newedge->c = (float) (dx * s1->coord.x + dy * s1->coord.y + (dx * dx + dy * dy) * 0.5);
  64.  
  65.   if (adx > ady) {              /* covers case dy = 0.0 */
  66.     newedge->a = (float) 1.0;
  67.     newedge->b = (float) (dy / dx);
  68.     newedge->c /= (float) dx;
  69.   }
  70.   else {                        /* covers case dx = 0.0 */
  71.     newedge->b = (float) 1.0;
  72.     newedge->a = (float) (dx / dy);
  73.     newedge->c /= (float) dy;
  74.   }
  75.  
  76.   newedge->edgenbr = nedges;
  77.   out_bisector (newedge, imgIO, value);
  78.   nedges += 1;
  79.  
  80.   return (newedge);
  81. }
  82.  
  83.  
  84. struct Site *
  85. intersect (el1, el2)
  86.      struct Halfedge *el1, *el2;
  87. {
  88.   struct Edge *e1, *e2, *e;
  89.   struct Halfedge *el;
  90.   float d, xint, yint;
  91.   int right_of_site;
  92.   struct Site *v;
  93.  
  94.   e1 = el1->ELedge;
  95.   e2 = el2->ELedge;
  96.   if ((e1 == (struct Edge *) NULL) || (e2 == (struct Edge *) NULL))
  97.     return ((struct Site *) NULL);
  98.   if (e1->reg[1] == e2->reg[1])
  99.     return ((struct Site *) NULL);
  100.  
  101.   d = e1->a * e2->b - e1->b * e2->a;
  102.   if ((-1.0e-10 < d) && (d < 1.0e-10))
  103.     return ((struct Site *) NULL);
  104.  
  105.   xint = (e1->c * e2->b - e2->c * e1->b) / d;
  106.   yint = (e2->c * e1->a - e1->c * e2->a) / d;
  107.  
  108.   if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
  109.       (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
  110.        e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
  111.     el = el1;
  112.     e = e1;
  113.   }
  114.   else {
  115.     el = el2;
  116.     e = e2;
  117.   }
  118.  
  119.   right_of_site = xint >= e->reg[1]->coord.x;  /* Boolean */
  120.   if ((right_of_site && el->ELpm == le) ||
  121.       (!right_of_site && el->ELpm == re))
  122.     return ((struct Site *) NULL);
  123.  
  124.   v = (struct Site *) getfree (&sfl);
  125.   v->refcnt = 0;
  126.   v->coord.x = xint;
  127.   v->coord.y = yint;
  128.  
  129.   return (v);
  130. }
  131.  
  132. /* returns 1 if p is to right of halfedge e */
  133. int
  134. right_of (el, p)
  135.      struct Halfedge *el;
  136.      struct Point *p;
  137. {
  138.   struct Edge *e;
  139.   struct Site *topsite;
  140.   int right_of_site, above, fast;
  141.   float dxp, dyp, dxs, t1, t2, t3, yl;
  142.  
  143.   e = el->ELedge;
  144.   topsite = e->reg[1];
  145.   right_of_site = p->x > topsite->coord.x;
  146.   if (right_of_site && (el->ELpm == le))
  147.     return (1);
  148.   if (!right_of_site && (el->ELpm == re))
  149.     return (0);
  150.  
  151.   if (e->a == 1.0) {
  152.     dyp = p->y - topsite->coord.y;
  153.     dxp = p->x - topsite->coord.x;
  154.     fast = 0;
  155.     if ((!right_of_site & e->b < 0.0) | (right_of_site & e->b >= 0.0)) {
  156.       above = dyp >= e->b * dxp;
  157.       fast = above;
  158.     }
  159.     else {
  160.       above = p->x + p->y * e->b > e->c;
  161.       if (e->b < 0.0)
  162.         above = !above;
  163.       if (!above)
  164.         fast = 1;
  165.     }
  166.  
  167.     if (!fast) {
  168.       dxs = topsite->coord.x - (e->reg[0])->coord.x;
  169.       above = e->b * (dxp * dxp - dyp * dyp) <
  170.         dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
  171.       if (e->b < 0.0)
  172.         above = !above;
  173.     }
  174.   }
  175.   else {                        /*e->b==1.0 */
  176.     yl = e->c - e->a * p->x;
  177.     t1 = p->y - yl;
  178.     t2 = p->x - topsite->coord.x;
  179.     t3 = yl - topsite->coord.y;
  180.     above = t1 * t1 > t2 * t2 + t3 * t3;
  181.   }
  182.  
  183.   return (el->ELpm == le ? above : !above);
  184. }
  185.  
  186.  
  187. void
  188. endpoint (e, lr, s, imgIO, value)
  189.      struct Edge *e;
  190.      int lr;
  191.      struct Site *s;
  192.      Image *imgIO;
  193.      int value;
  194. {
  195.   e->ep[lr] = s;
  196.   ref (s);
  197.  
  198.   if (e->ep[re - lr] == (struct Site *) NULL)
  199.     return;
  200.  
  201.   out_ep (e, imgIO, value);
  202.   deref (e->reg[le]);
  203.   deref (e->reg[re]);
  204.   makefree ((struct Freenode *) e, &efl);  /* diff. types, parm 1 !! */
  205. }
  206.  
  207.  
  208. float
  209. dist (s, t)
  210.      struct Site *s, *t;
  211. {
  212.   float dx, dy;
  213.  
  214.   dx = s->coord.x - t->coord.x;
  215.   dy = s->coord.y - t->coord.y;
  216.  
  217.   return ((float) sqrt (dx * dx + dy * dy));
  218. }
  219.  
  220.  
  221. void
  222. makevertex (v)                  /* S. Fortune orig. routine: type int */
  223.      struct Site *v;
  224. {
  225.   v->sitenbr = nvertices;
  226.   nvertices += 1;
  227.   out_vertex (v);
  228. }
  229.  
  230.  
  231. void
  232. deref (v)
  233.      struct Site *v;
  234. {
  235.   v->refcnt -= 1;
  236.   if (v->refcnt == 0)
  237.     makefree ((struct Freenode *) v, &sfl);
  238. }
  239.  
  240. void
  241. ref (v)
  242.      struct Site *v;
  243. {
  244.   v->refcnt += 1;
  245. }
  246.